Simplified engineering geomorphic unit-based seismic site characterization of the detailed area plan of Dhaka city, Bangladesh

Severe failure of improperly designed and poorly constructed structures may occur due to the amplified and prolonged ground motion during an earthquake, and hence prediction of the ground motion characteristics at the soil surface is crucial. In this study, based on the prepared simplified engineering geomorphic map, we performed a one-dimensional (1D) nonlinear site response analysis for seismic site characterization of the recently proposed Detailed Area Plan (DAP) area of Dhaka City, the Capital of Bangladesh. The engineering geomorphic unit-based map was prepared from image analysis and verified with the collected borehole data and surface geology map. The study area was classified into three major geomorphic units and seven sub-units subject to the subsurface soil profiles. Nine earthquake time histories, seven from the PEER NGA WEST2 data set and two synthetics, and seven identified subsurface soil profiles were used for nonlinear site response analysis, along with the BNBC 2020 uniform hazard spectrum as the target spectrum. For the selected earthquake ground motions, the near-surface soil response of the DAP area showed de-amplification of acceleration in the short period and amplification of acceleration in the long period. The amplified long-period acceleration could cause severe damage in inappropriately designed and poorly constructed long-period structures. The outcome of this study could be used to prepare a seismic risk-sensitive land use plan for the future development of the DAP of Dhaka City.

www.nature.com/scientificreports/ shows the earthquakes (with a magnitude of Mw ≥ 5) located in and around Bangladesh between 1908 and 2015.
The tectonic belt of the Himalayan system is expected to experience large earthquakes between M 7.5 and M 8.5 as the Indian plate moves northward at a rate of 4 cm/year and northeastward at a rate of 6 cm/year 45 . The Dauki Fault, which is thought to have caused the 1897 Great Indian Earthquake 46 , was activated three times in the last thousand years, according to recent paleo-seismological research [47][48][49] . As a result, the rate of plate movement and regular occurrence of large-magnitude earthquakes along this plate boundary suggests that Bangladesh, Nepal, Bhutan, Myanmar, and the northern and northeastern parts of India, are seismically active regions 43 . The Chittagong-Tripura fold belt has a noticeably different landscape due to the N-S trending anticlines, which resulted from the collision of the Indian and the Burmese plates, resting on the décollement megathrust of the arc-trench system of the Indian plate and the Burma Arc 50 . This zone is a consequence of the ongoing deformation along the subduction zone. In addition, the Neogene deformation system of the Chittagong-Tripura fold belt has minor seismicity indicating a state of locked condition of the plates 44 .

Method
This study prepared a simplified engineering geomorphic map of the Dhaka DAP area. Then, over two hundred boreholes with various lab test results and downhole seismic test data are collected to verify the simplified engineering geomorphic unit and delineate unit-wise engineering properties. Later, nine strong ground motion data are selected to execute site response analysis. Supplementary Fig. S2 shows the methodological framework adopted for this study. Further details of each component are explained in the subsequent sub-sections.
Simplified engineering geomorphic map. To reduce earthquake risk, an engineering geomorphic map is essential for zoning any area. It is a land classification map based on the geomorphic unit, which helps to know the soil type variation on the surface. We have used both boreholes' data and satellite images to prepare the engineering geomorphic map. The Landsat 5 image of 1988 and the Landsat 8 image for 2020 from USGS Earth Explorer are used for this research and are corrected using radiometric and atmospheric correction. Using the image from 1988, we digitized the basic geographic unit through visual interpretation (using natural and false colors). We then verified the geomorphic unit map from boreholes data and surface geology data from GSB (Geological Survey of Bangladesh). Later, NDWI (Normalized Difference Water Index) detected water bodies and landfills for both images of 1988 and 2020. We extract the water bodies from images by NDWI, give values 0 to 1, and then analyze the changes of water bodies in the images. From the studied images, landfill areas and water bodies are detected for the present time. Then these two units are added to the previously prepared map.
Engineering soil properties. Primary and secondary data have been collected to delineate subsurface engineering properties. Primary data includes Dhaka University Master Plan and Development of Open Space Management System to Response Scenario Earthquake in Dhaka Metropolitan Area projects data as some authors were directly involved in acquiring, processing, and interpreting Standard Penetration Test (SPT) borehole, geotechnical lab test, and downhole seismic test.
The secondary data are collected from the Geological Survey of Bangladesh (GSB), Creative Soil Investigations, Prosoil, Comprehensive Disaster Management Programme (CDMP). All data contain information about soil geotechnical properties, SPT N value, lithological properties, grain size, liquid limit, plastic limit, shear wave velocity (Vs), cohesion, friction angle, etc.
We have identified three major types of soil profiles from the collected borehole data considering the geomorphic units. They are Pleistocene terrace (overconsolidated clay or silty clay or clayey silt (OC)), Holocene alluvium (HA), and Landfill (LAN). Further, OC is divided into OC 1, OC 2, and OC 3 based on the subsurface soil layer. Similarly, HA and LAN are divided as HA 1, HA 2, and LAN-OC LAN-HA, respectively. LAN-OC and LAN-HA mean LAN underlay by OC and HA, respectively.
Site response analysis. Shear wave velocity of the study area. Shear wave velocity (Vs) of the near-surface soils is critical in seismic site response analysis 51,52 . Although it is costly and not simple like SPT, it has advantages over SPT-N value 53 . For example, measurements can be taken in hard soils where the SPT and Cone Penetration Test (CPT) are impossible or not permitted or in areas where undisturbed samples are difficult to collect. In addition, it is a fundamental mechanical property of soil. Therefore, it is directly related to the small-strain shear modulus, a required parameter in analytical procedures for estimating dynamic soil response, site response analysis, and soil-structure interaction analyses 54 . The study's Downhole Seismic Test (DST) locations are shown in Fig. 2. The DST is a direct method to estimate the Vs. and is considered the most reliable and accurate 55  where, h i and v i denote the thickness (m) and Vs (m/s) of the ith layer, respectively, and N denotes the number of layers in the top 30 m depth.
The shear wave velocity profiles of all geomorphic units are given in Supplementary Fig. S4.
(1) V S = 92.1 × N 0.337 (All soils) www.nature.com/scientificreports/ Shear modulus reduction and material damping curves. Shear modulus and soil material damping curves are prerequisites for nonlinear site response analysis 56 . As shear modulus and damping ratio reduction curves are not available for Dhaka city, it is a prevalent practice to use standard curves from the literature. The widely used curves for categorizing dynamic soil behavior are represented in many studies [57][58][59] . In this study, to fit the damping ratio and shear modulus reduction curve, as this empirical model is one of the most acceptable among practitioners for determining the variability of soft soil nonlinear properties, are used from Darendeli 57  yy : natural pressure. If the ratio becomes more than one, the soil is overconsolidated; equal to one usually means consolidated, and less than one implies under-consolidation 60 . The PI of soils was collected from laboratory tests. The number of loading cycles in the paper used is ten recommended by the DEEPSOIL manuals. The effective vertical stress of the soil layer was calculated by multiplying the soil layer thickness with the unit weight of the soil. If the soil is saturated, then pore water pressure will be minus from the vertical stress. Then, the effective vertical stress becomes σ = γ * h − γ w * h ; where h is thickness and γ: unit weight of soil, γ w : unit weight of water 61 . Using the soil PI, OCR, and the number of loading cycles, this model can predict the baseline of the damping ratios and the modulus reduction curves. Supplementary Fig. S5 illustrated the reference 57 , fit and current curves of shear modulus, damping ratio, and shear strength.
Strong ground motion selection. The strong ground motion time history (acceleration, velocity, and displacement) data is a prerequisite for dynamic response analysis of soil and structure. However, selecting strong ground motion is challenging as it depends not only on epicentral distance and rupture length but also on other factors, such as subsurface geology, directivity, etc. In the case of a moderate seismicity-prone and developing country, Bangladesh, the lack of strong ground motion data (Mw > 7.0 earthquake did not occur for more than 100 years in Bangladesh) is becoming more challenging. Therefore, we have used the most widely used spectral matching technique in this study to select the strong ground time history acceleration data. Firstly, we have downloaded strong motion raw data (horizontal component, H1) from the PEER NGA WEST2 database considering the similar response spectra to the Bangladesh National Building Code (BNBC) 2020 response spectra (target response spectra) for Maximum Credible Earthquake (MCE) and SB type soil condition (Vs 360 to 800 m/s; Supplementary Fig. S6) 36 . The details of the downloaded earthquake accelerograms are given in Supplementary Table S1. Afterward, we adjusted earthquake accelerograms by matching them with BNBC 2020 SB_MCE response spectra 36 using the wavelets algorithm proposed by Atik and Abrahamson 41 . Note that we have used SeismoMatch for spectral matching. Supplementary Fig. S7 shows the matched and target response spectra (BNBC 2020 SB_MCE). In addition, for this study, we have generated two artificial earthquake accelerograms matched to target response spectra (BNBC 2020 SB_MCE) by using SeismoArtif software. Both far-field inter and intra-plate earthquakes have been considered to generate synthetic accelerograms. In summary, seven earthquake accelerograms from seven earthquakes and two synthetic accelerograms have been considered for site response analysis. Supplementary Fig. S8 illustrates the detailed (step-by-step) procedure of ground motion selection for the Kobe 1995 earthquake.
Computing one-dimensional surface site response analysis. When motion is propagated vertically and horizontally, 1D site response analysis evaluates the soil column's impacts on earthquake time histories vertically. Frequency domain and time-domain analysis are the most often used techniques in 1D analysis. However, timedomain response analysis is more realistic than frequency-domain analysis as it consider changes in soil characteristics over time and the propagation of input motion into the soil 62 . The nonlinear time domain site response analysis was performed based on two model backbone curves and the hysteresis behavior of soil in terms of shear stress-strain loading-unloading.
The dynamic response of multi-degree-of-freedom systems exposed to base excitation is determined by nonlinear site response. This study utilizes DEEPSOILv7 software 23 for calculating nonlinear response analysis. Determining nonlinear, linear, and equivalent linear response analysis in DEEPSOIL is possible. Several soil constitutive models are accessible in DEEPSOIL, including the MKZ model with pressure-dependent behavior 63 and the new General Quadratic/Hyperbolic (GQ/H) model with non-masing requirements 64 . Numerous models for reference fit curves for shear modulus reduction and damping ratio curves are also available in this computer program. Among all the models, the GQ/H model is used as it is updated, widely used, and has the best fit for sandy and clayey soils. The GQ/H model for nonlinear site response analysis is based on the physics of wave propagation in heterogeneous media. The GQ/H model assumes that a quadratic-hyperbolic function can represent the stress-strain relationship of the soil. This quadric model used two lines of initial shear strength and shear strength at failures that capture small and large strain behavior of soil into a continuous curve because within linear boundaries initial shear strength and shear strength at failure intersect at some reference shear strain and create backbone curve 64,65 . The equation for the backbone curve of the GQ/H model is given below: www.nature.com/scientificreports/ τ is the shear stress, τ max is the shear strength at failure, γ is the shear strain, γ r is the reference shear strain and θ r is a curve fitting parameter used to adjust the curves of shear stress and strain curves and has no effect on boundary conditions. The shear strength at failure was measured using Mohr-Coulomb failure Criteria. Mohr-Coulomb failure criteria is given below Here, τ max is the shear strength of soil, c is the cohesion of soil, σ ′′ is effective stress, and φ is the frictional angle of soil. The cohesion and frictional angle are measured by lab test and given in Table 1. The effective stress measurement is already mentioned above. From Eq. (5), estimated shear strength was used as a user-defined parameter in DEEPSOIL for the GQ/H model 66 . The backbone curve uses Konder 67 proposed values as the reference shear strain 64,65 . The GQ/H model attempts to become flexible with minor strain soil behavior while simultaneously introducing shear strength failure. After the backbone curve created from the GQ/H model, the curve incorporates with hysteresis behavior of soil. In the hysteresis model, the shear stress-strain curve follows the backbone curve as τ = F bb (γ ) for initial shear loading. Therefore, the backbone curve of GQ/H model becomes The nonlinear time domain site response uses the extended unloading -loading masing rules to model hysteresis behavior. However, using masing rules in hysteresis damping calculation can overestimate damping in large strain. Therefore, Philips and Hashash 68 proposed a modulus reduction and damping factor (MRDF) approach. When this MRDF approach is used with the GQ/H model, the damping reduction matches well with the reference curves of Darendeli 65 . The MRDF reduction factor is where, p 1 , p 2 and p 3 are nondimensional parameters selected to obtain the best fit with the target damping curve. After defining the soil model parameters, the reduction factor parameters are chosen to best fit the damping curve. The computer program already has a built-in correlation of the GQ/H model and MRDF parameters for the best-fitted damping curve 65 . The MDRF model was developed for non-masing rules in the cyclic behavior of soil 69 . As observed in laboratory tests, the damping behavior of soils can be better captured by a non-Masing hysteresis model that shrinks the size of hysteresis loops. The used non-Masing hysteresis model considers that the soil's response may change over time as it undergoes deformation and stress. This model can capture the damping behavior more accurately by shrinking the size of the hysteresis loop, which represents the difference in energy input and output over one loading cycle.
Computing amplification factor (AF). Site coefficient or Amplification Factor (AF) is commonly used to express ground motion amplification at different spectral periods. We have used Eq. (8) to estimate the AF at different spectral periods.
where A is the amplification value of T second, R soil is the response of soil on the surface at T second, and R inputmotion is the input motion of earthquake at T second 70,71 .

Results
Engineering geomorphic map. The engineering geomorphic unit map of the DAP of the Dhaka City area was prepared and verified with the boreholes and surface geological data. The region has been divided into three simplified geomorphic units. They are the Holocene deposit (HA), Pleistocene deposit (OC), and Landfill (LAN) (Fig. 2). Nearly 50% of the area comprises Holocene deposits, and 33% includes Pleistocene deposits. The sticky, stiff to very stiff silty clay or clayey silt and reddish brown to yellowish brown color is the unique characteristics of OC underlain by medium dense to dense silty sand 34 . HA is mainly composed of light grey to dark grey color, very soft to soft clayey silt or silty clay, and gray to brown, loose to dense sandy silt or silty sand 34 . The common mineral components of the HA unit are quartz, feldspar, and mica 35,37,72 . The LAN unit is characterized by grey loose silty sand to sand and soft clayey silt. River dredging sand-filling sites are very vulnerable to liquefaction hazards; for example, the 1995 Kobe earthquake caused massive liquefaction in Port Island's reclaimed soil deposits. www.nature.com/scientificreports/ formation (Fig. 3). The identified two prominent verities of HA are HA 1 and HA 2. HA 1 consists of, from top to bottom, Holocene clay, Holocene silt, and Holocene sand. However, HA 2 consists of HA 1 and OC 1, where OC 1 is overlain by HA 1 (Fig. 3). In addition, there are two types of LAN soil units where LAN overlies OC (LAN OC) and LAN underlain by HA (LAN HA). The average filling depth in these two types is 3-4 m, and most of the LAN was composed of loose, unconsolidated river sand. Landfill sand and Holocene sand are susceptible to liquefaction. The liquefaction phenomenon will appear on the surface or not; it depends on the relative thickness of the liquefiable and non-liquefiable layer 73 . In addition, the Liquefaction Potential Index (LPI) value is also a good and widely used criterion for estimating surface manifestations of liquefaction hazards 74 . Engineering properties of the identified seven soil profiles of three geomorphic units are given in Table 1.
Site response analysis. The geomorphic unit-based nonlinear surface response spectra for nine selected input ground motions are illustrated in Fig. 4. The response spectrums indicate various outcomes depending on the geomorphic unit and input earthquakes ground motion; however, the overall observation is that spectral acceleration is de-amplified (ground motion amplitude decreases compared to input ground motion) in the short period and amplified (amplitude of ground motion increase compared to input ground motion) in the long period. All the responses are measured considering a 30 m soil column where Vs is 360-450 m/s at the bottom of the soil column. As we know, earthquake ground motion on the ground surface depends on the source, path, and soil characteristics of selected sites. The ground motion with definite frequency interacts with the predominant frequency of soil; as a result, a variation in the response spectrum on the surface is observed. Figure 4, represents the response spectral acceleration of OC, HA, and LAN units. The OC unit has three subunits namely: OC 1, OC 2, and OC3. Figure 4a-c represent the response of OC 1, OC 2, and OC3, respectively. The OC soil (OC 1, OC 2, and OC 3) (Fig. 4a-c) represents that the responses are de-amplified in a short period (0.1 s) while peak acceleration at around 0.7 s which means the response is amplified in the long period. However, similar responses are observed from OC 1, OC 2, and OC 3, as the Madhupur Clay and Silt engineering properties are close (Table 1). Figure 4d,e show the response spectrum for LAN-HA and LAN-OC, respectively. As expected, responses from LAN units (LAN-HA and LAN-OC) are not identical. The LAN-HA shows predominant peak acceleration at 0.2 s (Fig. 4d), whereas the LAN-OC unit show peak response at around 0.5 s (Fig. 4e). The upper 3 m of the LAN unit is filled with loose river sand (Table 1). After the 3 m depth, the LAN-OC unit has Madupur clay, silt, and Dupitila layer and LAN-HA has Holocene clay, silt, sand, and Dupitila formation. This means that two LAN units distinguish soil types and properties. These different soil properties and layer thicknesses cause different responses for LAN units. The last unit is the HA unit, HA 1 and HA 2. HA 1 means the soil up to 30 m entirely consists of Holocene clay, silt, and sand, whereas HA 2 consists of both Holocene clay, silt, sand, Madupur clay, silt, and Dupitila sand. The spectral acceleration of HA 1 and HA 2 illustrated in Fig. 4f,g, respectively. Figure 4f,g show slightly different response spectra than other units. The predominant period of the considered soil column ranges from approximately 0.40 s to 0.50 s, and the input ground motion predominant period is mainly concentrated around 0.40 s (from 0.22 to 0.48 s), expected Loma Prieta (around 0.22 s). Therefore, although Loma Prieta has the highest matched PGA (0.31 g), its response is less. On the other hand, Kobe (peak period 0.46 s) has matched PGA of 0.27 g and shows a higher response than Loma Prieta. The predominant period of OC 1, OC 2, OC 3, and LAN-OC is around 0.40 s; therefore, their responses are similar as expected. On the other hand, the predominant period of HA 1, HA 2, and LAN-HA is around 0.50 s, showing similar responses. This analysis calculates response spectra concerning frequency-independent viscous damping ratio using MRDF curve fitting for Darendali's curve fitting factor 68 . Figure 4h illustrates the AF of seven geomorphic subunits at different spectral periods for Gorkha earthquake input. The AF is decreased (< 1) in short period (before approximately 0.5 s) and increased (> 1) in a long period (after approximately 0.5 s). www.nature.com/scientificreports/ The layer-wise nonlinear response for OC 1 unit is illustrated in Fig. 5. As expected, the top layer, Madhupur clay, has a maximum response spectrum than Madhupur silt and Dupitila Formation (Fig. 5c) compared with the input earthquake motion (Gorkha 2015), the input motion has maximum spectral acceleration at 0.17 s with 0.66 g. In contrast, the responses have maximum spectral acceleration at 0.53 s with 0.83 g (Fig. 5c). The responses are de-amplified in a short period until 0.1 s. Figure 5d shows the difference between the linear and nonlinear responses of Madhupur clay. Spectral acceleration from the linear analysis is overestimated in the short period and underestimated in the long period. However, the nonlinear response shows de-amplification in a short period and amplification in a long period. As most of the soil of the DAP area is relatively soft, the response values would be overestimated using linear analysis because the linear model cannot effectively predict ground motion during the short period when peak shear strain is between 0.01 and 0.1% 13 .

Discussion
In this study, site responses are estimated for seven types of engineering geomorphic sub-units under three main units using nonlinear analysis. We observed that the response spectrum shifted towards the long period with maximum acceleration more significant than the input ground motion (Figs. 4 and 5). The soft/loose sandy (LAN (mainly loose sand) (SPT N value 0-10)) and clayey (HA (SPT N value 0-10) and Madhupur clay (SPT N value 5-25)) soil up to 30 m depth is responsible for amplifying acceleration over a long period and de-amplification of acceleration in short periods. However, according to the National Earthquake Hazard Reduction Program    www.nature.com/scientificreports/ (NEHRP) site AF, for 30 m soil depth, the response is amplified in a short period 52 . On the other hand, the soft/loose sandy and clayey soil of DAP of Dhaka city area does not follow the NHERP site AF and shows deamplification over a short period (Fig. 4h).
In a nonlinear analysis, the higher frequencies are amplified at high shear moduli, whereas the lower frequencies are amplified at low shear moduli (Supplementary Fig. S5). The de-amplification and amplification on the response spectrum depend on the shear wave velocity and the other soil properties. Hashash and Park 75 showed that because of the decreased shear wave velocity at the top 70 m of the profile, the HA profile exhibits less amplified acceleration in short periods than the Pleistocene profile and a slight shifting to longer periods.
De-amplification or amplification of the ground surface response spectrum mainly depends on input ground motion characteristics and soil properties (loose or dense, soft or stiff). The soil shows a response according to the frequency and amplitude of the ground motion. As we know, soft soil damps high-frequency ground motion and amplifies low-frequency ground motion. In this study, due to the presence of nonlinear soft clay, silt, and very loose sand deposits, the short-period ground motion de-amplified, and the long-period ground motion amplified. When the input motion frequency matches the soil's natural frequency, it amplifies the input motion.
The HA covers most of the study area. OC and an artificial LAN cover the rest of the site. The response spectrums are different according to soil type and shear wave velocity. We can see that OC with maximum acceleration over a long period than HA and LAN deposits. On the contrary, HA exhibits long-lasting amplification than OC. The artificial LAN shows a mixture of them. High frequencies are attenuated more than low frequencies when passing through a column of soft, loose soil that consists primarily of sandy soil. High frequencies are amplified more in OC sediments than in thick and soft HA due to the higher natural frequencies of OC sediments. Thicker Holocene sediments have lower natural frequencies, resulting in amplification in the lower frequencies www.nature.com/scientificreports/ than the OC sediments. Dhaka City is very susceptible to far-field earthquakes. The low-frequency tremor with adequate amplitude can cause significant damage to low-frequency structures due to resonance effects. In summary, for the considered earthquakes, the near-surface soil response of DAP of Dhaka city shows de-amplification of acceleration in the short period and amplified acceleration in the long period. This type of spectrum is controlled by the characteristics and properties of soft clayey and sandy soil, and response spectra differ according to the geomorphic unit. For example, the Holocene soil unit accelerates toward a long period more than the Pleistocene unit. However, the Pleistocene deposits exhibit higher acceleration values than those of the Holocene.

Conclusion
This paper examined the engineering geomorphic unit-based nonlinear site response of the DAP area of Dhaka City. First, the simplified engineering geomorphic unit map was prepared by image analysis and then verified using borehole data and the surface geology map of GSB. The identified three major geomorphic units are OC, HA, and LAN. Further, OC, HA, and LAN are divided into OC 1, OC 2, OC 3, HA 1, HA 2, LAN-OC, and LAN-HA sub-units considering the subsurface soil layer.
In total, nine input ground motions (seven from PEER NGA WEST2 data set and two synthetics) were used for 1D nonlinear site response analysis of seven geomorphic sub-units of DAP of Dhaka city. The overall observation from the nonlinear soil site response analysis is that the responses are de-amplified in the short period while amplified in the long period.
OC's spectral acceleration amplitude is higher than HA units. However, HA units show longer peak spectral acceleration duration than OC units. This phenomenon might cause the resonance of tall buildings on HA units, consequently increasing the probability of damage. As expected, LAN-OC units show similar behavior of OC and LAN-HA to HA. Dhaka city's vigorous vertical and lateral expansion has caused it to create a large built environment over soft soil. Therefore, the findings of this study can be used as an essential indicator to prepare a seismic risk-sensitive land use plan for the future development of the DAP of Dhaka City to reduce the seismic risk. A similar procedure could be adopted for other earthquake-vulnerable cities to understand the seismic risk.

Data availability
Strong Ground Motion Data are available at https:// ngawe st2. berke ley. edu/. The SPT and DST data from mentioned organizations could be appropriately acquired by contact with them.